#include <windows.h>
#include "config.h"

#define PI 3.14159265358979323846264338327950288419716939937510

#define a 124367
#define c 57634
#define m 2147483648
//unsigned int x=1234;

__device__ float rnd2(unsigned int idx)
{
    unsigned int x=1234;
    int aa=a;
    int res = 1;
    int t = idx+1;
    while (t)
    {
		if (t & 1)
			res *= aa;
		t >>= 1;
    }
    x = (res*x)%m;
	return ((float) x)/m;
}

__device__ float rnd(unsigned int &_x)
{
    //unsigned int x=1234;
	//for (int i = 0; i < skip; i++)    
    //    _x=(a*_x+c)%m;
	
	_x=(a*_x+c)%m;
	return ((float) _x)/m;
	//return ((float)rand())/RAND_MAX;
}

__device__ float f_device(float x, float y)
{
	return (exp( sin(PI * x) * cos(PI * y) ) + 1.0) / 256.0;
}

__global__ void MCCVKernel(int x1, int y1, int x2, int y2, int chunkSize, float *S)
{

	//long i = blockIdx.x * blockDim.x + threadIdx.x;
	//long j = blockIdx.y * blockDim.y + threadIdx.y;
	
	long ind = blockIdx.x * blockDim.x + threadIdx.x;
	
	//long id = i*N + j;		
	//S[id] = 0;
	
	//float *valV = new float[chunkSize], *valW = new float[chunkSize];
	__shared__ float valV[CHUNKSIZE];
	__shared__ float valW[CHUNKSIZE];
	//float valV[CHUNKSIZE];
	//float valW[CHUNKSIZE];

	float denom = (x2-x1)*(y2-y1);
	float rnd_x = 0, rnd_y = 0;//, rnd_z = 0;
	float Fxy = 0;//, Fxy_max = 0;
	float result = 0.0f;
	float corVW=0.0f, meanV=0.0;
	
	//just W(X1, X2) = (X1 + X2), Mean(W) = (X1max-X1min)/2 + (X2max-X2min)/2
	//float meanW = (x2-x1)/2.0f + (y2-y1)/2.0f, sigmaW2=0.0f;
	float meanW = 0.5f*(x1+x2+y2+y1), sigmaW2=0.0f;
	
	unsigned int num = ind*chunkSize*2;
		
	for (int i = 0; i < chunkSize; i++)
	{	
		rnd_x = rnd2(num+2*i) * (x2-x1)+x1;
		rnd_y = rnd2(num+2*i+1) * (y2-y1)+y1;
		Fxy = f_device(rnd_x, rnd_y);

		valV[i] = Fxy;
		valW[i] = rnd_x + rnd_y;	
		////meanWExper += valW[i];

		sigmaW2 += (valW[i] - meanW)*(valW[i] - meanW);
		meanV += valV[i];		
	}	
	
	
	meanV /= chunkSize;
	for (int i = 0; i < chunkSize; i++)
		corVW += (valV[i] - meanV) * (valW[i] - meanW);
		
	
	float alpha = corVW / sigmaW2;
	
	
	for (int i = 0; i < chunkSize; i++)
	{
		result += (valV[i] - alpha*(valW[i] - meanW)) * denom;
	}
	

	result /= chunkSize;
	
	S[ind] = result;
}

int MCCVCuda(int x1, int y1, int x2, int y2, int nsamples, int chunkSize, float *output, LARGE_INTEGER &time)
{
	//dim3 threadsPerBlock(16, 16);
	//dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
	float *S;
	
	int chunksCount = nsamples / chunkSize;
	
	cudaMalloc(&S, chunksCount * sizeof(float));	
	//MCCVKernel<<<chunksCount, chunkSize>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<1, chunkSize>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<nsamples/15, 15>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<chunksCount/15, 15>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<2, 15>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<chunksCount/60, 60>>>(x1, y1, x2, y2, chunkSize, S);
	MCCVKernel<<<8, chunksCount/8>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<1, chunksCount>>>(x1, y1, x2, y2, chunkSize, S);
	//MCCVKernel<<<chunkSize, 1>>>(x1, y1, x2, y2, chunkSize, S);

    cudaThreadSynchronize();
    
    QueryPerformanceCounter(&time);
	cudaMemcpy(output, S, chunksCount*sizeof(float), cudaMemcpyDeviceToHost);
	cudaFree(S);
	return 0;
} 